! THIS CODE IS ONLY USED FOR THE ALTERNATIVE MODEL WITH WITHIN-OCCUPATION COMPLEMENTARITIES (REPORTED IN APPX F.5)
MODULE PARAMS
implicit none
!integer, parameter :: n_a=100,n_z=20,n_x=30,n_id=1000,n_time=400,maxiterx=200,load=1,keeplaborfixed=0,maxiterk=200
integer, parameter :: n_a=100,n_z=21,n_tau=5,n_x=14,n_id=1000,n_time=400,maxiterx=20,load=1,keeplaborfixed=0&
,maxiterk=1,maxiterxcalib=1,n_iterminw=1
integer, parameter :: n_d=n_z*n_tau
double precision, parameter :: critdiff=1E-7,critxcalib=1E-3
double precision alpha,beta,eta,r,delta,rho(2),gamma,labor(n_d,2),pitot,empfirm(n_d),distfirm(n_z,n_tau),&
avect(n_a),zvect(n_z),tauvect(n_tau),xvect(n_x),cmin,sigma,rhoext,pselect,wage(n_d,n_x)&
,unemprate(n_d,n_x),minw(n_x),skewz,skewx,kurtx,b_unemp,bmat(n_d,n_x),pd(n_d,n_d),theta
integer firmconvert(n_z,n_tau),zcomp(n_d),taucomp(n_d)

double precision , dimension(n_a,n_d,n_x):: distw,f1,vinit
double precision , dimension(n_a,n_x):: coef,vinside,vsum
double precision :: laborfixed(n_d,2),mufixed(n_d,n_x)


double precision , dimension(n_a,n_d,n_x,2):: vlast 
integer , dimension(n_a,n_d,n_x,2):: iprmat

double precision mean_z,z_inf,z_sup,updown,delta_z,std_z,z_left,z_right,rhoz,std_y
double precision mean_x,x_inf,x_sup,delta_x,std_x,x_left,x_right,rhox,std_tempx



CONTAINS

subroutine dss_search(index_i,index_z,index_h,index_u,vmax,imax,lower,upper)
implicit none
integer index_i,index_z,index_h,index_u,lower,upper,imax,ii,iitemp
double precision vmax,maxtemp

  imax=lower
  maxtemp=-1E+20

  do ii=lower,upper
  if (vfun(ii,index_i,index_z,index_h,index_u)>maxtemp) then
  iitemp=ii
  maxtemp=vfun(ii,index_i,index_z,index_h,index_u)
  endif
  
  enddo
  imax=iitemp
  vmax=maxtemp
end subroutine


function vfun(i_choice,i_state,index_z,index_x,index_u)
implicit none
integer i_choice,i_state,index_z,index_x,z_choice,index_u
double precision cons,vfun,util,fnext,inc


if (index_u==1) then
inc=bmat(index_z,index_x)
else if (index_u==2) then
inc=wage(index_z,index_x)
endif


cons=inc+avect(i_state)*(1.+r)-avect(i_choice)+pitot



if (cons>cmin) then
util=cons**(1.-gamma)/(1.-gamma)
else
util=cmin**(1.-gamma)/(1.-gamma)-1000.d0*abs(cmin-cons)**2.d0
endif


fnext=f1(i_choice,index_z,index_x)


vfun=util+beta*fnext   
                                     
end function

double precision function choiceprob(i_state,index_z,index_x)
implicit none

integer i_state,index_z,index_x
double precision nom,denom


choiceprob=exp(vinit(i_state,index_z,index_x)/(sigma*rhoext)-coef(i_state,index_x))/vinside(i_state,index_x)

if (vinside(i_state,index_x)==0.d0) then

if  (vinit(i_state,index_z,index_x)/(sigma*rhoext)==coef(i_state,index_x)) then
choiceprob=1.d0
else
choiceprob=0.d0
endif
endif

if (choiceprob>1.d0) then
print *, 'tfun>1',i_state,index_z,index_x,choiceprob,vinit(i_state,index_z,index_x),&
vinit(i_state,index_z,index_x)/(sigma*rhoext),coef(i_state,index_x),vinside(i_state,index_x)
pause
endif


if (choiceprob.ne.choiceprob) then
print *, 'in choiceprob', nom,denom,vinit(i_state,index_z,index_x),vinside(i_state,index_x)
pause
endif


end function





double precision function normalcdf(x)
implicit none
double precision x,inside,erf,a1,a2,a3,a4,a5,rat,p

p = 0.3275911d0
a1 = 0.254829592d0
a2 = -0.284496736d0
a3 = 1.421413741d0
a4 = -1.453152027d0
a5 = 1.061405429d0

inside=(abs(x)-0.d0)/(sqrt(2.d0))
rat=1.d0/(1.d0+p*inside)

erf=1-(a1*rat+a2*rat**2+a3*rat**3+a4*rat**4+a5*rat**5)*exp(-inside**2)
if (x<0) erf=-erf

normalcdf=0.5d0*(1+erf)

end function

      
ENDMODULE PARAMS
      
program solvemodel
use params
implicit none
      
double precision laborfirm(n_d,2),kfirm(n_d),yfirm(n_d),labornew(n_d,2),avgl(n_d),sdl(n_d),mufirm(n_d,n_x),&
sdyfirm,sdempfirm,pz(n_z,n_z),px(n_x),&
distznew(n_z),diffdistz,ksup,laborold(n_d,2),minwclose(2),unempall(2),sdwage,kdem,wtarget(n_x),minwavgwagecollar(2)&
,disttau(n_tau),distz(n_z)
double precision tau_inf,tau_sup,mean_tau,std_tau,delta_tau,tau_left,tau_right

!taumin,ptail
double precision, dimension(n_d,n_x):: wagenew,mudem,mudemnew,musup,wageold

double precision, dimension(n_a,n_d,n_x)::distwnew,distwselect,distwnoselect,f1new,changedecrate,changefirmrate

!px(n_x,n_x),std_tempx,rhox,distgvnew(n_z,n_x),distgvold(n_z,n_x)


double precision sum0,sum1,sum2,sum3,sum4,sum5,grad,difflabor(2),diffdist,diffa,temp,maxtemp,flow,nothing(4),rate,ptemp,&
mufirmtemp,laborfirmtemp(2),diffmu,z,x,diffwage(2),autocorrw,autocovw,diffdistrel,difflaborrel,diffwagerel,diffmurel,decilez(n_z),&
diffk,minwvect(5),avgwage(n_x),sum6,sum7,sum8,sum9,sum10,autocovy,autocorry,autocovemp,autocorremp,diffxcalib(n_x),&
unempxspec(n_x),minwclosexspec(n_x),ptemp0,avgwage_level(n_x),diffminw(2)

integer i,j,h,hh,jj,ii,iterdist,iterz,stopz,iterx,stay,itermu,stopmu,n,k,kk,iord,stopdist,id,t,iterk,minwgp(5),iterxcalib,&
jz,jtau,jjz,jjtau,l,iterminw

double precision, dimension(n_a,n_d,n_x):: choiceprobmat

 double precision zbar,zgrid,ppz,qqz
 double precision, dimension(n_z,n_z,n_z)::ptempz1,ptempz2,ptempz3,ptempz4,panz,pantempz
 
 double precision xbar,xgrid,ppx,qqx
 double precision, dimension(n_x,n_x,n_x)::ptempx1,ptempx2,ptempx3,ptempx4,panx,pantempx
 
 integer ord(n_a),invord(n_a),mini,maxi,ilo(n_a),ihi(n_a),stoplo,stophi,glo,ghi
 
 integer, dimension(n_id,n_time):: a_sim,z_sim,x_sim,u_sim,choose_sim,tau_sim,d_sim
 double precision, dimension(n_id,n_time):: val_sim,inc_sim,c_sim
 double precision, dimension(n_id,n_time):: randomz,randomx,randomk,randomu,randomchoose
      

 
!double precision, dimension(:,:,:,:,:), ALLOCATABLE :: choiceprobmat
!double precision, dimension(:,:,:,:), ALLOCATABLE :: choiceprobforcemat
!integer AllocateStatus, DeAllocateStatus

open(unit=903,file="load_pars.txt", action='read')
    READ(903,*) gamma
    READ(903,*) beta
    READ(903,*) alpha
    READ(903,*) eta
    READ(903,*) delta   
    READ(903,*) rho(1)
    READ(903,*) rho(2) 
    READ(903,*) skewz 
    READ(903,*) std_z
    READ(903,*) rhoz  
!    READ(903,*) skewx
!    READ(903,*) kurtx
!    READ(903,*) mean_x
!    READ(903,*) std_x
!    READ(903,*) rhox    
    READ(903,*) pselect
    READ(903,*) sigma
    READ(903,*) rhoext
    READ(903,*) mean_tau
    READ(903,*) std_tau
close(903)

open(unit=903,file="load_theta.txt", action='read')
    READ(903,*) theta
close(903)


open(unit=2010,file="r.txt")
read (2010,*) r
close(2010)


open(unit=90,file="minw.txt", action='read')
    READ(90,*) minwvect(1)
    READ(90,*) minwvect(2)
    READ(90,*) minwvect(3)
    READ(90,*) minwvect(4)
    READ(90,*) minwvect(5)
close(90)

!minwvect(2)=minwvect(1)*60.08d0/50.67d0

!minwvect(2)=minwvect(1)*1.05

ii=int(n_x/5)
write(*,*) ii
write(*,*) minwvect

minwgp=0;minw=0.d0

do i=1,n_x
!k=1
!do while (i>ii*k.and.k<=5)
!k=k+1
!enddo
!minwgp(i)=k
!minw(i)=minwvect(k)
if (i<=7) then
minwgp(i)=1
minw(i)=minwvect(1)
else
minwgp(i)=2
minw(i)=minwvect(2)
endif

enddo



!do i=1,n_x
!minwgp(i)=floor(i*1.d0/5)+1 
!minw(i)=minwvect(minwgp(i))
!enddo

do i=1,n_x
write(*,*) minw(i),minwgp(i)!,minwvect(minwgp(i))
enddo
!pause


b_unemp=minwvect(3)

!in the counterfactual, maintain the same b as in the bm
!b_unemp=40.d0/55.d0*420.


      open(unit=27,file="b_unemp.txt")
      write (27,*) b_unemp
      close(27)

!sigma=0.01d0
!rhoext=1.d0

cmin=0.00001d0


!!!! When we want to control for the unemployment affect in welfare analysis,
!!!! we will just give them the wage they would earn without the minw.    

bmat=b_unemp
!open(unit=909,file="wage_nominw.txt")
!    do j=1,n_z
!    read(909,*) (bmat(j,h),h=1,n_x)
!    enddo
!close(909)



!!! Start the guesses for labor stats !!!
labor=1.5d0
pitot=0.d0



!minw=1.d0
wage=1.d0
unemprate=0.d0
f1=0.d0



avect=0.d0
if (n_a>1) then
! avect(n_a)=10000.d0
 avect(n_a)=100.d0
 
      open(unit=87,file="ahigh.txt")
      read (87,*) avect(n_a)
      close(87)



 do i=2,n_a
 avect(i)=avect(1)+(i-1.)*(avect(n_a)-avect(1))/(n_a-1)
 enddo
 

 
! do i=1,n_a
! temp=.2d0
! avect(i)=1.04d0**(i-1)-1.d0
! enddo

endif
!write(*,*) avect
!pause

!!! Divide and conquer algorithm (Gordon and Qiu) for the asset choice !!!
ord=0
ord(1)=1;ord(n_a)=2
invord(1)=1;invord(2)=n_a


i=3
do while (minval(ord)==0)

mini=1
do while (ord(mini+1)>0)
mini=mini+1
enddo


maxi=mini
do while (ord(maxi)==0.or.mini.eq.maxi)
maxi=maxi+1
enddo
if (maxi>mini+1) then
ord(int((mini+maxi)/2))=i
invord(i)=int((mini+maxi)/2)
endif
i=i+1
!write(*,*) ord
!write(*,*) mini,maxi
!pause
enddo

ilo=0
ihi=0

ilo(n_a)=1

do i=2,n_a-1
ii=i;stoplo=0
do while (ii>=1.and.stoplo==0)
if (ord(ii)<ord(i)) then
ilo(i)=ii
stoplo=1
else
ii=ii-1
endif
enddo

ii=i;stophi=0
do while (ii<=n_a.and.stophi==0)
if (ord(ii)<ord(i)) then
ihi(i)=ii
stophi=1
else
ii=ii+1
endif
enddo

enddo

!!!!!!!! z distribution of firms !!!!!!
mean_z=-10.d0

if (rho(1)>0.2d0) mean_z=0.d0
!if (rho(1)>0.2d0) mean_z=0.d0

updown=3.
std_y = std_z/ SQRT(1 - rhoz**2)
z_inf = mean_z - updown*std_y
z_sup = mean_z + updown*std_y
delta_z = updown*std_y / (n_z - 1d+0)


do i=1,n_z
   zvect(i) = exp(z_inf + (z_sup - z_inf) * (i-1) / (n_z - 1))
end do

   do i=1,n_z
   z_left  = (log(zvect(1)) + delta_z - rhoz*log(zvect(i))-(1.-rhoz)*mean_z)/ std_z
   z_right = (log(zvect(n_z))- delta_z - rhoz*log(zvect(i))-(1.-rhoz)*mean_z)/ std_z
   pz(i,1) = normalcdf(z_left)
   pz(i,n_z) = 1d+0 - normalcdf(z_right)
   do j=2,n_z-1
     z_right = (log(zvect(j)) + delta_z - rhoz*log(zvect(i))-(1.-rhoz)*mean_z)/ std_z
     z_left  = (log(zvect(j)) - delta_z - rhoz*log(zvect(i))-(1.-rhoz)*mean_z)/ std_z
     pz(i,j) = (normalcdf(z_right) - normalcdf(z_left) ) !/prob_man_y
   end do
   end do
   
  !!!! ROUWENHORST !!!!!
!  
!         
!
!      zbar=(((n_z-1)/(1-rhoz**2.d0))**0.5d0)*std_z
!     
!      zvect(1)=mean_z-zbar
!      zgrid=zbar/((n_z-1.d0)/2.d0)
!      do j=2,n_z
!      zvect(j)=zvect(j-1)+zgrid
!      enddo
!      
!      !!! Skew is a parameter that dictates the thickness of the right tail in the ergo dist. 
!      !!! Has to be in (-(1-rhoz)/(1+rhoz),(1-rhoz)/(1+rhoz))
!      
!      ppz=(rhoz+1)/2*(1.+skewz)!+0.1d0
!      qqz=(rhoz+1)/2*(1.-skewz)!-0.1d0
!       
!      panz(1,1,1)=1.d0
!
!      do n=2,n_z
!
!      ptempz1=0.d0
!      ptempz2=0.d0
!      ptempz3=0.d0
!      ptempz4=0.d0
!
!      do i=1,n-1
!      do j=1,n-1
!      ptempz1(i,j,n)=panz(i,j,n-1)
!      ptempz2(i,j+1,n)=panz(i,j,n-1)
!      ptempz3(i+1,j,n)=panz(i,j,n-1)
!      ptempz4(i+1,j+1,n)=panz(i,j,n-1)
!      enddo
!      enddo
!      pantempz(:,:,n)=ppz*ptempz1(:,:,n)+(1-ppz)*ptempz2(:,:,n)+(1-qqz)*ptempz3(:,:,n)+qqz*ptempz4(:,:,n)
!
!      panz(1,:,n)=pantempz(1,:,n)
!      panz(n,:,n)=pantempz(n,:,n)
!      do i=2,n-1
!      panz(i,:,n)=pantempz(i,:,n)/2
!      enddo
!
!      enddo
!      pz=panz(:,:,n_z)      
!      do j=1,n_z
!      zvect(j)=exp(zvect(j))
!      enddo


!disttau=1./n_tau 
!tauvect(1)=0.2d0;tauvect(2)=0.1d0;tauvect(3)=-0.1d0;tauvect(4)=-0.20d0
!taumin=1.d0
!ptail=2.d0 

  
!taumin=0.42d0
!temp=0.90d0
!
!tauvect(n_tau)=taumin*(1.-temp)**(-1./ptail)
!
!rate=0.75d0
!
!ptemp0=temp/((1.-rate**n_tau)/(1.-rate))
!
!
!tauvect(1)=taumin*(1.-ptemp0)**(-1./ptail)
!do j=2,n_tau
!tauvect(j)=taumin*((taumin/tauvect(j-1))**ptail-ptemp0*rate**(j-1))**(-1./ptail)
!enddo
!
!disttau(1)=1.-(taumin/tauvect(1))**ptail 
!do j=2,n_tau
!disttau(j)=(taumin/tauvect(j-1))**ptail-(taumin/tauvect(j))**ptail
!enddo
!
!tauvect=tauvect-1.d0  

updown=2.d0
tau_inf = mean_tau - updown*std_tau
tau_sup = mean_tau + updown*std_tau
delta_tau = updown*std_tau / (n_tau - 1d+0)

do i=1,n_tau
   tauvect(i) = exp(tau_inf + (tau_sup - tau_inf) * (i-1) / (n_tau - 1))
end do

   tau_left  = (log(tauvect(1)) + delta_tau-mean_tau)/ std_tau
   tau_right = (log(tauvect(n_tau))- delta_tau-mean_tau)/ std_tau
   disttau(1) = normalcdf(tau_left)
   disttau(n_tau) = 1. - normalcdf(tau_right)
   do j=2,n_tau-1
     tau_right = (log(tauvect(j)) + delta_tau - mean_tau)/ std_tau
     tau_left  = (log(tauvect(j)) - delta_tau - mean_tau)/ std_tau
     disttau(j) = normalcdf(tau_right) - normalcdf(tau_left)
   end do


sum0=sum(disttau)
disttau=disttau/sum0
tauvect=tauvect-1.


 
   
distz=1.d0/n_z

iterz=1;stopz=0
do while (iterz<2000.and.stopz==0)

do jj=1,n_z
sum0=0.d0
do j=1,n_z
sum0=sum0+pz(j,jj)*distz(j)
enddo
distznew(jj)=sum0
enddo

diffdistz=maxval(abs(distz-distznew))
distz=distznew
if (diffdistz<0.0000001d0) then
stopz=1
else
iterz=iterz+1
endif
enddo


do l=1,n_tau
do j=1,n_z
distfirm(j,l)=distz(j)*disttau(l)
enddo
enddo

!write(*,*) distz
!pause


!write(*,*) sum(distz),sum(disttau),sum(distfirm)
!pause


!if (load==1) then
!      open(unit=292,file="distfirm.txt")
!      do j=1,n_z
!      read (292,*) (distfirm(j,l),l=1,n_tau)
!      enddo
!      close(292)
!else
!distfirm=0.d0
!distfirm(1,:)=1.d0
!endif

i=1
do l=1,n_tau
do jz=1,n_z
firmconvert(jz,l)=i
zcomp(i)=jz
taucomp(i)=l
i=i+1
enddo
enddo



do j=1,n_d
do jj=1,n_d
if (taucomp(j)==taucomp(jj)) then
pd(j,jj)=pz(zcomp(j),zcomp(jj))
else
pd(j,jj)=0.d0
endif
enddo
enddo


      open(unit=2945,file="tauvect.txt")
      do i=1,n_tau
      write (2945,*) tauvect(i)
      enddo
      close(2945)

      open(unit=2955,file="disttau.txt")
      do i=1,n_tau
      write (2955,*) disttau(i)
      enddo
      close(2955)

      open(unit=294,file="zvect.txt")
      do i=1,n_z
      write (294,*) zvect(i)
      enddo
      close(294)
                 
      open(unit=297,file="pz.txt")
      do j=1,n_z
      write (297,*) (pz(j,jj),jj=1,n_z)
      enddo
      close(297)
      
      open(unit=292,file="distfirm.txt")
      do j=1,n_z
      write (292,*) (distfirm(j,l),l=1,n_tau)
      enddo
      close(292)
      

do j=1,n_z
if (j<=6) then
decilez(j)=1
else if (j<=8) then
decilez(j)=2
else if (j<=9) then
decilez(j)=3
else if (j<=10) then
decilez(j)=4
else if (j<=11) then
decilez(j)=6
else if (j<=12) then
decilez(j)=7
else if (j<=13) then
decilez(j)=8
else if (j<=15) then
decilez(j)=9
else
decilez(j)=10
endif
enddo



!!! x process for workers !!!
!updown=kurtx*2.
!
!std_tempx = std_x/ SQRT(1 - rhox**2)
!x_inf = mean_x - updown*std_tempx
!x_sup = mean_x + updown*std_tempx
!delta_x = updown*std_tempx / (n_x - 1d+0)
!
!
!do i=1,n_x
!   xvect(i) = exp(x_inf + (x_sup - x_inf) * (i-1) / (n_x - 1))
!end do
!
!   do i=1,n_x
!   x_left  = (log(xvect(1)) + delta_x - rhox*log(xvect(i))-(1.-rhox)*mean_x)/ std_x
!   x_right = (log(xvect(n_x))- delta_x - rhox*log(xvect(i))-(1.-rhox)*mean_x)/ std_x
!   px(i,1) = normalcdf(x_left)
!   px(i,n_x) = 1d+0 - normalcdf(x_right)
!   do j=2,n_x-1
!     x_right = (log(xvect(j)) + delta_x - rhox*log(xvect(i))-(1.-rhox)*mean_x)/ std_x
!     x_left  = (log(xvect(j)) - delta_x - rhox*log(xvect(i))-(1.-rhox)*mean_x)/ std_x
!     px(i,j) = (normalcdf(x_right) - normalcdf(x_left) ) !/prob_man_y
!   end do
!   end do
!   
!
!
!
!do j=1,n_x
!
!do jj=min(j+1,n_x),n_x
!px(j,jj)=px(j,jj)*(1.-skewx)
!enddo
!
!px(j,1)=px(j,2);px(j,n_x)=px(j,n_x-1)
!
!enddo
!
!do j=1,n_x
!px(j,:)=px(j,:)/sum(px(j,:))
!enddo
!
!write(*,*) (sum(px(j,:)),j=1,n_x)

!!!! Blue collar !!!
!px(1)=0.14;px(2)=0.15;px(3)=0.14;px(4)=0.14;px(5)=0.13;px(6)=0.11;px(7)=0.09;px(8)=0.06;px(9)=0.03;px(10)=0.01
!!!! White collar !!!
!px(11)=0.04;px(12)=0.03;px(13)=0.03;px(14)=0.05;px(15)=0.07;px(16)=0.09;px(17)=0.13;px(18)=0.18;px(19)=0.22;px(20)=0.17
!
!wtarget(1)=67.d0;wtarget(2)=71.d0;wtarget(3)=74.d0;wtarget(4)=78.d0;wtarget(5)=81.d0;wtarget(6)=84.d0;wtarget(7)=88.d0;&
!wtarget(8)=95.d0;wtarget(9)=108.d0;wtarget(10)=143.d0
!
!wtarget(11)=74.d0;wtarget(12)=80.d0;wtarget(13)=84.d0;wtarget(14)=88.d0;wtarget(15)=91.d0;wtarget(16)=97.d0;wtarget(17)=104.d0;&
!wtarget(18)=114.d0;wtarget(19)=128.d0;wtarget(20)=157.d0
!
!sum0=sum(px(1:10));sum1=wtarget(1)
!do j=1,10
!px(j)=px(j)/sum0*0.63d0
!wtarget(j)=wtarget(j)/sum1
!enddo
!sum0=sum(px(11:n_x))
!do j=11,n_x
!px(j)=px(j)/sum0*0.37d0
!wtarget(j)=wtarget(j)/sum1
!enddo

!!! Blue collar !!!
px(1)=0.098d0;px(2)=0.105d0;px(3)=0.098d0;px(4)=0.098d0;px(5)=0.091d0;px(6)=0.077d0;px(7)=0.063d0



!!! White collar !!!
px(8)=0.0206d0;px(9)=0.0288d0;px(10)=0.037d0;px(11)=0.053d0;px(12)=0.074d0;px(13)=0.090d0;px(14)=0.070d0


wtarget(1)=67.d0;wtarget(2)=71.d0;wtarget(3)=74.d0;wtarget(4)=78.d0;wtarget(5)=81.d0;wtarget(6)=84.d0;wtarget(7)=88.d0

wtarget(8)=88.d0;wtarget(9)=91.d0;wtarget(10)=97.d0;wtarget(11)=104.d0;wtarget(12)=114.d0;wtarget(13)=128.d0;wtarget(14)=157.d0

sum0=sum(px)
px=px/sum0

sum1=wtarget(1)
wtarget=wtarget/sum1

!!!! If we want to match the minW baseline levels of wages for each x !!!!
!wtarget(1)=1.d0;wtarget(2)=1.0594793d0;wtarget(3)=1.104167;wtarget(4)=1.163787;wtarget(5)=1.208522;wtarget(6)= 1.253287d0
!wtarget(7)=1.312950d0
!wtarget(8)=1.310331d0;wtarget(9)=1.357892d0;wtarget(10)=1.447257d0;wtarget(11)=1.551644d0;wtarget(12)=1.700704d0;wtarget(13)=1.909405d0
!wtarget(14)=2.342760d0

xvect=1.d0
!xvect(1)=1.d0
!do j=2,10
!xvect(j)=xvect(j-1)*1.d0
!enddo
!
!xvect(11)=1.d0
!do j=12,n_x
!xvect(j)=xvect(j-1)*1.0d0
!enddo


distw=0.d0
do h=1,n_x
distw(1,:,h)=px(h)/(1.d0*n_d)
enddo


!!!! ROUWENHORST !!!!!

!      xbar=kurtx*(((n_x-1)/(1-rhox**2.d0))**0.5d0)*std_x
!     
!      xvect(1)=mean_x-xbar
!      xgrid=xbar/((n_x-1.d0)/2.d0)
!      do j=2,n_x
!      xvect(j)=xvect(j-1)+xgrid
!      enddo
!      
!      !!! Skew is a parameter that dictates the thickness of the right tail in the ergo dist. 
!      !!! Has to be in (-(1-rhoz)/(1+rhoz),(1-rhoz)/(1+rhoz))
!      
!      ppx=(rhox+1)/2*(1.+skewx)!+0.1d0
!      qqx=(rhox+1)/2*(1.-skewx)!-0.1d0
!       
!      panx(1,1,1)=1.d0
!
! 
!      do n=2,n_x
!
!      ptempx1=0.d0
!      ptempx2=0.d0
!      ptempx3=0.d0
!      ptempx4=0.d0
!
!      do i=1,n-1
!      do j=1,n-1
!      ptempx1(i,j,n)=panx(i,j,n-1)
!      ptempx2(i,j+1,n)=panx(i,j,n-1)
!      ptempx3(i+1,j,n)=panx(i,j,n-1)
!      ptempx4(i+1,j+1,n)=panx(i,j,n-1)
!      enddo
!      enddo
!      pantempx(:,:,n)=ppx*ptempx1(:,:,n)+(1-ppx)*ptempx2(:,:,n)+(1-qqx)*ptempx3(:,:,n)+qqx*ptempx4(:,:,n)
!
!      panx(1,:,n)=pantempx(1,:,n)
!      panx(n,:,n)=pantempx(n,:,n)
!      do i=2,n-1
!      panx(i,:,n)=pantempx(i,:,n)/2
!      enddo
!
!      enddo
!      px=panx(:,:,n_x) 
!
!     
!      do j=1,n_x
!      xvect(j)=exp(xvect(j))
!      enddo

!      open(unit=294,file="xvect.txt")
!      do i=1,n_x
!      write (294,*) xvect(i)
!      enddo
!      close(294)
      
                 
      open(unit=292,file="px.txt")
      do j=1,n_x
      write (292,*) px(j)
      enddo
      close(292)




if (load==1.or.load==2) then


open(unit=905,file="labor.txt", action='read')
    do j=1,n_d
    READ(905,*) labor(j,1),labor(j,2)
    enddo
close(905)



open(unit=907,file="pi.txt", action='read')
    READ(907,*) pitot
close(907)




open(unit=909,file="wage.txt")
open(unit=910,file="unemp.txt")
    do j=1,n_d
    read(909,*) (wage(j,h),h=1,n_x)
        read(910,*) (unemprate(j,h),h=1,n_x)
    enddo
close(909)
close(910)



open(unit=9112,file="xvect.txt")
    do h=1,n_x
    read(9112,*) xvect(h)
    enddo
close(9112)


if (load==1) then

open(unit=911,file="f1.txt")
    do i=1,n_a
    do j=1,n_d
    read(911,*) (f1(i,j,h),h=1,n_x)
    enddo
    enddo
close(911)


open(unit=908,file="distw.txt")
    do i=1,n_a
    do j=1,n_d
    read(908,*) (distw(i,j,h),h=1,n_x)
    enddo
    enddo
close(908)

endif

!!!!When recalibrating x's to match the same levels as in the baseline !!!!
!xvect(1)=xvect(1)*0.999


do j=1,n_d
do h=1,n_x
musup(j,h)=sum(distw(:,j,h))
enddo
enddo

endif


!Start dens from 0!
!    do j=1,n_d
!    do h=1,n_x
!    i=1
!    sum0=sum(distw(:,j,h))
!    distw(i,j,h)=sum0
!    do i=2,n_a
!    distw(i,j,h)=0.d0
!    enddo
!    enddo
!    enddo



!xvect=xvect*0.998d0


mufixed=musup
laborfixed=labor

write(*,*) 'starting'

iterminw=1;diffminw=100.d0
do while (iterminw<=n_iterminw)

iterxcalib=1;diffxcalib=100.d0
do while (iterxcalib<=maxiterxcalib)

ksup=100.d0;diffk=100.d0;iterk=1
do while (iterk<=maxiterk)

diffwagerel=100.d0;difflaborrel=100.d0;diffa=100.d0;iterx=0
do while ((difflaborrel>critdiff.or.diffdistrel>critdiff).and.iterx<maxiterx)



glo=1;ghi=n_a
 !$OMP PARALLEL private(j,h,k,iord,i,glo,ghi)
 !$OMP DO SCHEDULE(RUNTIME) ORDERED 
do j=1,n_d
do h=1,n_x
do k=1,2
  do iord=1,n_a
    i=invord(iord)
    
  if (ilo(i)==0.or.iprmat(ilo(i),j,h,k)==n_a) then
  glo=1
  else
  glo=iprmat(ilo(i),j,h,k)
  endif  
  
  if (ihi(i)==0) then
  ghi=n_a
  else
  ghi=iprmat(ihi(i),j,h,k)
  endif 
  
  call dss_search(i,j,h,k,vlast(i,j,h,k),iprmat(i,j,h,k),glo,ghi)

! call dss_search(i,j,h,k,vlast(i,j,h,k),iprmat(i,j,h,k),1,n_a)

 
 enddo
 enddo  
 enddo
 enddo  
 !$OMP END DO
 !$OMP END PARALLEL
 
! j=7;h=2;k=2
! do i=1,n_a
!! write(*,*) zcomp(j),i,iprmat(i,j,h,k),iprmat(ilo(i),j,h,k),iprmat(ihi(i),j,h,k),avect(iprmat(i,j,h,k))/avect(i)-1.
! write(*,*) zcomp(j),i,iprmat(i,j,h,k),avect(iprmat(i,j,h,k))/avect(i)-1.
!
! enddo
! pause
 

 !$OMP PARALLEL private(i,j,h)
 !$OMP DO SCHEDULE(RUNTIME) ORDERED 
do i=1,n_a
do j=1,n_d
do h=1,n_x

vinit(i,j,h)=unemprate(j,h)*vlast(i,j,h,1)+(1.-unemprate(j,h))*vlast(i,j,h,2)

enddo  
enddo
enddo  
 !$OMP END DO
 !$OMP END PARALLEL



  !$OMP  PARALLEL DO PRIVATE(i,j,h,maxtemp)

do i=1,n_a
do h=1,n_x

maxtemp=-1E+20 !.d0
!meantemp=0.d0

  !$OMP  PARALLEL DO PRIVATE(j) SHARED(i,h) REDUCTION(max : maxtemp)
do j=1,n_d
maxtemp=max(vinit(i,j,h),maxtemp)
enddo

  !$OMP END PARALLEL DO
 coef(i,h)=maxtemp/(sigma*rhoext)
 

enddo
enddo
!$OMP END PARALLEL DO



vinside=0.d0;vsum=0.d0
   !$OMP  PARALLEL DO SCHEDULE(RUNTIME)  PRIVATE(i,j,h,sum0,flow)

do i=1,n_a
do h=1,n_x

sum0=0.d0
  !$OMP  PARALLEL DO &
  !$OMP  PRIVATE(j,flow) SHARED(i,h)      &
  !$OMP  REDUCTION(+ : sum0)

do j=1,n_d
flow=exp(vinit(i,j,h)/(sigma*rhoext)-coef(i,h))
sum0=sum0+flow
enddo

  !$OMP END PARALLEL DO
vinside(i,h)=sum0

vsum(i,h)=sigma*log(sum0**rhoext)+sigma*rhoext*coef(i,h)

enddo
enddo
  !$OMP END PARALLEL DO
  

  
  
!$OMP  PARALLEL DO SCHEDULE(RUNTIME)  PRIVATE(i,j,h,sum0,jj)
do i=1,n_a
do j=1,n_d
do h=1,n_x

sum0=0.d0
  !$OMP  PARALLEL DO &
  !$OMP  PRIVATE(jj) SHARED(i,j,h)      &
  !$OMP  REDUCTION(+ : sum0) 
do jj=1,n_d
sum0=sum0+pd(j,jj)*vinit(i,jj,h)
enddo
  !$OMP END PARALLEL DO
f1new(i,j,h)=(1.-pselect)*sum0+pselect*vsum(i,h)
enddo
enddo
enddo

  !$OMP END PARALLEL DO
  
 !$OMP  PARALLEL DO SCHEDULE(RUNTIME)  PRIVATE(i,j,h)

do i=1,n_a
do j=1,n_d
do h=1,n_x
choiceprobmat(i,j,h)=choiceprob(i,j,h)

!if (choiceprobmat(i,j,h)>2.d0.or.choiceprobmat(i,j,h)<-0.1) then
!write(*,*) 'choiceprobmat',choiceprobmat(i,j,h),&
!vinit(i,j,h)/(sigma*rhoext),coef(i,h),vinside(i,h)
!pause
!endif

enddo
enddo
enddo
  !$OMP END PARALLEL DO

!write(*,*) 'choiceprobmat2',minval(choiceprobmat),maxval(choiceprobmat),sum(choiceprobmat)/n_a/n_x
!pause

!stopdist=0;iterdist=0
!do while (stopdist==0)
 

!distwselect=0.d0!;distwnoselect=0.d0
! !$OMP  PARALLEL DO PRIVATE(i,j,h,ii,jj,hh,temp,k)
!do jj=1,n_z
!do hh=1,n_x
! !$OMP  PARALLEL DO SHARED(jj,hh) PRIVATE(i,j,h,ii,temp,k)
!do i=1,n_a 
!do j=1,n_z 
!do h=1,n_x
!do k=1,2
!ii=iprmat(i,jj,h,k)
!temp=unemprate(jj,h)
!if (k==2) temp=1.-unemprate(jj,h)
!distwnoselect(ii,jj,hh)=distwnoselect(ii,jj,hh)+distw(i,jj,h)*temp*px(h,hh)*(1.-pselect)
!enddo
!enddo
!enddo
!  !$OMP END PARALLEL DO
!enddo
!enddo
!  !$OMP END PARALLEL DO
  
 
 distwnew=0.d0 
!$OMP  PARALLEL DO PRIVATE(i,j,ii,jj,hh,temp,k)
do jj=1,n_d
do hh=1,n_x
!$OMP  PARALLEL DO SHARED(jj,hh) PRIVATE(i,j,ii,temp,k)

do j=1,n_d
do i=1,n_a 
do k=1,2
ii=iprmat(i,j,hh,k)
temp=unemprate(j,hh)
if (k==2) temp=1.-unemprate(j,hh)
distwnew(ii,jj,hh)=distwnew(ii,jj,hh)+&
         distw(i,j,hh)*temp*(pselect*choiceprobmat(ii,jj,hh)+(1.-pselect)*pd(j,jj))

enddo
enddo
enddo
!$OMP END PARALLEL DO
enddo
enddo
!$OMP END PARALLEL DO
   

! distwnew=0.d0 
!!$OMP  PARALLEL DO PRIVATE(i,j,ii,jj,hh,temp,k,sum0)
!do ii=1,n_a
!do jj=1,n_d
!do hh=1,n_x
!sum0=0.d0
!  !$OMP  PARALLEL DO &
!  !$OMP  PRIVATE(i,j,k,temp) SHARED(ii,jj,hh)      &
!  !$OMP  REDUCTION(+ : sum0) 
!do j=1,n_d
!do i=1,n_a 
!do k=1,2
!if (ii==iprmat(i,j,hh,k)) then
!temp=unemprate(j,hh)
!if (k==2) temp=1.-unemprate(j,hh)
!sum0=sum0+distw(i,j,hh)*temp*(pselect*choiceprobmat(ii,jj,hh)+(1.-pselect)*pd(j,jj))
!endif
!enddo
!enddo
!enddo
!!$OMP END PARALLEL DO
!distwnew(ii,jj,hh)=sum0
!
!enddo
!enddo
!enddo
!!$OMP END PARALLEL DO
!   
   
! !$OMP  PARALLEL DO PRIVATE(i,j,h)
!  do i=1,n_a
!  do j=1,n_z
!  do h=1,n_x
!  distwnew(i,j,h)=distwselect(i,j,h)+distwnoselect(i,j,h)
!  enddo
!  enddo
!  enddo
!  !$OMP END PARALLEL DO
diffdistrel=0.d0
do i=1,n_a 
do j=1,n_d
do h=1,n_x
if (diffdistrel<abs((distw(i,j,h)-distwnew(i,j,h))/((distw(i,j,h)+distwnew(i,j,h))/2.d0+0.001))) then
diffdistrel=max(diffdistrel,abs((distw(i,j,h)-distwnew(i,j,h))/((distw(i,j,h)+distwnew(i,j,h))/2.d0+0.001)))
endif
enddo
enddo
enddo 


diffdist=sum(abs(distw-distwnew))





!if (diffdist<0.000001d0.or.iterdist>0) then
!stopdist=1
!else
!distw=distwnew
!iterdist=iterdist+1
!endif
!
!enddo !stopdist!
 
!write(*,*) 'dist iter', iterdist,diffdist,stopdist
!do j=1,n_z
!write(*,*) (distw(1,j,h),h=1,n_x)
!write(*,*) (distwnew(1,j,h),h=1,n_x) 
!enddo
!pause



if (abs(sum(distwnew)-1.d0)>0.000001d0.or.maxval(distwnew)>1.d0.or.minval(distwnew)<0.d0) then
write(*,*) 'distw problem',sum(distw),sum(distwnew),maxval(distwnew),minval(distwnew),sum(choiceprobmat)/n_a/n_x
pause
endif
grad=0.95d0
distw=(1.-grad)*distwnew+grad*distw
f1=f1new


do j=1,n_d
do h=1,n_x
musup(j,h)=sum(distw(:,j,h))
enddo
enddo

laborold=labor
wageold=wage


mudem=musup
labor=0.d0
do j=1,n_d
sum0=0.d0
sum1=0.d0
do h=1,n_x
if (minwgp(h)==1) then
sum0=sum0+mudem(j,h)**rho(1)*xvect(h)**rho(1)
else 
sum1=sum1+mudem(j,h)**rho(2)*xvect(h)**rho(2)
endif

enddo
labor(j,1)=sum0**(1./rho(1))
labor(j,2)=sum1**(1./rho(2))

enddo

if (keeplaborfixed==1) then
labor=laborfixed
mudem=mufixed
endif


do j=1,n_d
jz=zcomp(j);jtau=taucomp(j)
laborfirm(j,:)=labor(j,:)/distfirm(jz,jtau)
mufirm(j,:)=mudem(j,:)/distfirm(jz,jtau)
kfirm(j)=(alpha*eta/((1.+tauvect(jtau))*(r+delta))*zvect(jz)*laborfirm(j,1)**(theta*eta)*&
                    laborfirm(j,2)**((1.-alpha-theta)*eta))**(1./(1.-alpha*eta))
yfirm(j)=zvect(jz)*(kfirm(j)**alpha*laborfirm(j,1)**theta*laborfirm(j,2)**(1.-alpha-theta))**eta
empfirm(j)=sum(distw(:,j,:))/distfirm(jz,jtau)

enddo

pitot=0.d0
do j=1,n_d
jz=zcomp(j);jtau=taucomp(j)
pitot=pitot+(1.-eta)*yfirm(j)*distfirm(jz,jtau)
enddo
!pitot=0.d0 !!! CAREFUL; THIS IS JUST TEMPORARY

diffmu=1000.d0;itermu=0;stopmu=0
do while (stopmu==0)

do j=1,n_d
jz=zcomp(j);jtau=taucomp(j)
do h=1,n_x
x=xvect(h)
z=zvect(jz)/(1.+tauvect(jtau))
laborfirmtemp(1)=max(labor(j,1)/distfirm(jz,jtau),0.00001d0)
laborfirmtemp(2)=max(labor(j,2)/distfirm(jz,jtau),0.00001d0)

mufirmtemp=max(mudem(j,h)/distfirm(jz,jtau),0.00001d0)

if (minwgp(h)==1) then !bc!
wagenew(j,h)=theta*eta*(alpha*eta/(r+delta))**(alpha*eta/(1.-alpha*eta))*&
    z**(1./(1.-alpha*eta))*laborfirmtemp(1)**(theta*eta/(1.-alpha*eta)-rho(1))*&
   laborfirmtemp(2)**((1.-alpha-theta)*eta/(1.-alpha*eta))*x**rho(1)*mufirmtemp**(rho(1)-1.d0)

!write(*,*) 'in wagenew',z,x,laborfirmtemp,labor(j),distfirm(j),mufirmtemp
!pause  
    
if (wagenew(j,h)<minw(h)) then
!! Updated CC !!
if (rho(1)<0.98d0) then
mudemnew(j,h)=(theta*eta*(alpha*eta/(r+delta))**(alpha*eta/(1.-alpha*eta))*&
    z**(1./(1.-alpha*eta))*laborfirmtemp(1)**(theta*eta/(1.-alpha*eta)-rho(1))*&
    laborfirmtemp(2)**((1.-alpha-theta)*eta/(1.-alpha*eta))*x**rho(1)/minw(h))**(1./(1.-rho(1)))*distfirm(jz,jtau)
else
mudemnew(j,h)=0.d0
endif
    
else
mudemnew(j,h)=mudem(j,h)
endif

else !wc!
wagenew(j,h)=(1.-theta)*eta*(alpha*eta/(r+delta))**(alpha*eta/(1.-alpha*eta))*&
    z**(1./(1.-alpha*eta))*laborfirmtemp(1)**(theta*eta/(1.-alpha*eta))*&
   laborfirmtemp(2)**((1.-alpha-theta)*eta/(1.-alpha*eta)-rho(2))*x**rho(2)*mufirmtemp**(rho(2)-1.d0)

!write(*,*) 'in wagenew',z,x,laborfirmtemp,labor(j),distfirm(j),mufirmtemp
!pause  
    
if (wagenew(j,h)<minw(h)) then
!! Updated CC !!
if (rho(2)<0.98d0) then
mudemnew(j,h)=((1.-theta)*eta*(alpha*eta/(r+delta))**(alpha*eta/(1.-alpha*eta))*&
    z**(1./(1.-alpha*eta))*laborfirmtemp(1)**(theta*eta/(1.-alpha*eta))*&
    laborfirmtemp(2)**((1.-alpha-theta)*eta/(1.-alpha*eta)-rho(2))*x**rho(2)/minw(h))**(1./(1.-rho(2)))*distfirm(jz,jtau)
else
mudemnew(j,h)=0.d0
endif

else
mudemnew(j,h)=mudem(j,h)
endif

endif !bc wc!




enddo
enddo



labornew=0.d0
do j=1,n_d
sum0=0.d0
sum1=0.d0
do h=1,n_x
!! Updated CC !!
if (minwgp(h)==1) then
sum0=sum0+mudemnew(j,h)**rho(1)*xvect(h)**rho(1)
else
sum1=sum1+mudemnew(j,h)**rho(2)*xvect(h)**rho(2)
endif

enddo
labornew(j,1)=sum0**(1./rho(1))
labornew(j,2)=sum1**(1./rho(2))

enddo

diffmurel=0.d0
do j=1,n_d 
do h=1,n_x
if (diffmurel<abs((mudem(j,h)-mudemnew(j,h))/((mudem(j,h)+mudemnew(j,h))/2.d0+0.001))) then
diffmurel=max(diffmurel,abs((mudem(j,h)-mudemnew(j,h))/((mudem(j,h)+mudemnew(j,h))/2.d0+0.001)))
endif
enddo
enddo 

diffmu=maxval(abs(mudem-mudemnew))+maxval(abs(labor-labornew))
if (itermu==300.or.diffmurel<1E-10) then
stopmu=1
else
itermu=itermu+1

grad=0.!95d0

labor=grad*labor+(1.-grad)*labornew
mudem=grad*mudem+(1.-grad)*mudemnew
endif

!write(*,*) diffmu,itermu
!write(*,*) wagenew(3,:)
!pause 

!do j=1,n_d
!!write(*,*) labor(j,1),laborold(j,1),labor(j,2),laborold(j,2)
!write(*,*) itermu,wagenew(j,8),'lab',labor(j,2),mudemnew(j,h)
!
!enddo
!pause


enddo !!! Iter for mudem,labor





!j=3
!write(*,*) labor(j)
!write(*,*) (mudem(j,h),h=1,n_x)
!write(*,*) (wagenew(j,h),h=1,n_x)
!pause

wage=wagenew

do j=1,n_d
do h=1,n_x
unemprate(j,h)=(musup(j,h)-mudem(j,h))/musup(j,h)
if ((unemprate(j,h)<0.d0.or.unemprate(j,h)>1.d0).and.keeplaborfixed==0) then
write(*,*) 'unemp wrong',musup(j,h),mudem(j,h),wage(j,h),minw(h)
pause
endif

enddo
enddo

if (keeplaborfixed==1) then
unemprate=0.d0


if (maxval(abs(labor-laborfixed))>0.00000001d0.or.maxval(abs(mudem-mufixed))>0.0000001d0) then
pause 'not fixed'
endif
endif




sum0=0.d0
do i=1,n_a
do j=1,n_d
do h=1,n_x
sum0=sum0+distw(i,j,h)*avect(i)
enddo
enddo
enddo

diffa=abs((sum0-ksup)/ksup)
ksup=sum0

do kk=1,2
difflabor(kk)=maxval(abs(labor(:,kk)-laborold(:,kk)))/(sum(laborold(:,kk))/n_d)
diffwage(kk)=maxval(abs(wage(:,kk)-wageold(:,kk)))

enddo

!do j=1,10
!!write(*,*) labor(j,1),laborold(j,1),labor(j,2),laborold(j,2)
!write(*,*) wage(j,8),wageold(j,8),'lab',labor(j,2),laborold(j,2)!,'U',unemprate(j,h),musup(j,h),mudem(j,h)
!
!enddo
!pause


difflaborrel=0.d0
do kk=1,2
do j=1,n_d 
if (difflaborrel<abs((labor(j,kk)-laborold(j,kk))/((labor(j,kk)+laborold(j,kk))/2.d0+0.001))) then
difflaborrel=max(difflaborrel,abs((labor(j,kk)-laborold(j,kk))/((labor(j,kk)+laborold(j,kk))/2.d0+0.001)))
endif
enddo 
enddo

diffwagerel=0.d0
do j=1,n_d
do h=1,n_x
if (diffwagerel<abs((wage(j,h)-wageold(j,h))/((wage(j,h)+wageold(j,h))/2.d0+0.001))) then
!write(*,*) wage(j,h),wageold(j,h),j,h
diffwagerel=max(diffwagerel,abs((wage(j,h)-wageold(j,h))/((wage(j,h)+wageold(j,h))/2.d0+0.001)))
endif
enddo
enddo 

!pause


grad=0.5
!if (difflabor<0.1d0) grad=0.95d0

labor=grad*laborold+(1.-grad)*labor
wage=grad*wageold+(1.-grad)*wage


iterx=iterx+1


!write(*,*) 'avgxrho',avgxrho
!write(*,*) 'avgxrhonew',avgxrhonew
!write(*,*) 'labor',labornew
!write(*,*) (sum(distw(:,j,:)),j=1,n_z)
if (iterx-int(iterx/20)*20==0) then
!write(*,*) 'iterlabor',iterx,ksup,difflabor,'wages',diffwage,diffdistrel,diffmurel
!write(*,*) 'iterlabor',iterx,ksup,difflaborrel,'w',diffwagerel,diffdistrel,diffmurel
write(*,*) 'iterlabor',iterx,difflaborrel,'w',diffwagerel


!write(*,*) 'iterlabor 2',difflabor

endif
!pause

!pause

!endif

enddo


sum0=0.d0
do j=1,n_d
jz=zcomp(j);jtau=taucomp(j)
sum0=sum0+(1.-eta)*alpha*yfirm(j)/((1.+tauvect(j))*(r+delta))*distfirm(jz,jtau)
enddo
kdem=sum0

diffk=ksup/kdem-1.d0


if (abs(diffk)>0.001d0.and.iterk<maxiterk) then

r=r-0.01d0*diffk

if (r>0.1d0) then
r=0.1d0;iterk=5000
else if (r<0.001d0) then
r=0.001d0;iterk=6000
else
iterk=iterk+1
endif

else
iterk=maxiterk+100
endif

open(unit=2009,file="keqn.txt")
write (2009,*) ksup,kdem,r,iterk,diffk
close(2009)

open(unit=2010,file="r.txt")
write (2010,*) r
close(2010)

write(*,*) '---------',iterk,ksup,kdem,r,'------------'

enddo

do h=1,n_x
sum0=0.d0
do i=1,n_a
do j=1,n_d
sum0=sum0+distw(i,j,h)*wage(j,h)
enddo
enddo
avgwage(h)=sum0/px(h)
enddo

minwavgwagecollar=0.d0;sum0=0.d0
do i=1,7
sum0=sum0+avgwage(i)*px(i)
enddo
minwavgwagecollar(1)=minw(1)/(sum0/sum(px(1:7)))

sum0=0.d0
do i=8,n_x
sum0=sum0+avgwage(i)*px(i)
enddo
minwavgwagecollar(2)=minw(8)/(sum0/sum(px(8:n_x)))

avgwage_level=avgwage


sum1=avgwage(1)
avgwage=avgwage/sum1

diffxcalib=0.d0
if (iterxcalib<maxiterxcalib) then
diffxcalib(1)=0.d0
do h=2,n_x

diffxcalib(h)=avgwage(h)/wtarget(h)-1.d0


if (abs(diffxcalib(h))>critxcalib/n_x) then

if (h.ne.8) then

if (rho(minwgp(h))<0.d0) then 
xvect(h)=xvect(h)*(1.+0.5d0*diffxcalib(h))
else
xvect(h)=xvect(h)*(1.-0.5d0*diffxcalib(h))
endif

else

theta=theta*(1.+0.5d0*diffxcalib(h))

endif

endif

enddo

xvect=max(0.000000000001d0,xvect)
theta=max(0.01d0,min(1.-alpha-0.04,theta))

endif

if (sum(abs(diffxcalib))<critxcalib) then
iterxcalib=1000
else
iterxcalib=iterxcalib+1
endif

write(*,*) 'iterxcalib', iterxcalib, sum(abs(diffxcalib)),sum(abs(diffxcalib(1:7))),sum(abs(diffxcalib(9:n_x))),diffxcalib(8)

open(unit=304,file="avgwage.txt")
    do h=1,n_x
    write(304,*) xvect(h),avgwage(h),wtarget(h),avgwage_level(h)
    enddo
close(304)


open(unit=903,file="load_theta.txt")
    write(903,*) theta
close(903)


!write(*,*) '--------',iterxcalib
!do h=1,n_x,3
!write(*,*) xvect(h),avgwage(h),wtarget(h),diffxcalib(h)
!enddo
!pause

enddo

!j=1;
!write(*,*) 'mus'
!write(*,*) (mudem(j+1,h)/distfirm(j+1)/(mudem(j,h)/distfirm(j))-1,h=1,n_x,3)
!write(*,*) 'lambdas'
!write(*,*) (mudem(j+1,h)/mudem(j,h)-1,h=1,n_x,3)
!write(*,*) 'distws'
!write(*,*) (sum(distw(:,j+1,h))/sum(distw(:,j,h))-1,h=1,n_x,3)
!write(*,*) 'distw within z'
!write(*,*) (sum(distw(:,j+1,h))/sum(distw(:,j+1,:))/(sum(distw(:,j,h))/sum(distw(:,j,:)))-1,h=1,n_x,3)
!write(*,*) 'wages'
!write(*,*) (wage(j+1,h)/wage(j,h)-1,h=1,n_x,3)
!pause


!!! Stats for each firm of a given productivity !!!


sum0=0.d0;sum1=0.d0
do j=1,n_d
jz=zcomp(j);jtau=taucomp(j)
sum0=sum0+log(yfirm(j))**2.d0*distfirm(jz,jtau)
sum1=sum1+log(yfirm(j))*distfirm(jz,jtau)
enddo
sdyfirm=(sum0-sum1**2.d0)**0.5d0

sum0=0.d0;sum1=0.d0
do i=1,n_a
do j=1,n_d
do h=1,n_x
sum0=sum0+log(wage(j,h))**2.d0*distw(i,j,h)*(1.-unemprate(j,h))
sum1=sum1+log(wage(j,h))*distw(i,j,h)*(1.-unemprate(j,h))
enddo
enddo
enddo
sdwage=(sum0-sum1**2.d0)**0.5d0

sum0=0.d0;sum1=0.d0
do j=1,n_d
jz=zcomp(j);jtau=taucomp(j)
sum0=sum0+log(empfirm(j))**2.d0*distfirm(jz,jtau)
sum1=sum1+log(empfirm(j))*distfirm(jz,jtau)
enddo
sdempfirm=(sum0-sum1**2.d0)**0.5d0

sum0=0.d0;sum1=0.d0
do i=1,n_a
do j=1,n_d
do h=1,7
if (wage(j,h)<minw(h)*1.10d0) then
sum0=sum0+distw(i,j,h)*(1.-unemprate(j,h))
endif
sum1=sum1+distw(i,j,h)*(1.-unemprate(j,h))
enddo
enddo
enddo
minwclose(1)=sum0/sum1
unempall(1)=(sum(px(1:7))-sum1)/sum(px(1:7))

sum0=0.d0;sum1=0.d0
do i=1,n_a
do j=1,n_d
do h=8,n_x
if (wage(j,h)<minw(h)*1.10d0) then
sum0=sum0+distw(i,j,h)*(1.-unemprate(j,h))
endif
sum1=sum1+distw(i,j,h)*(1.-unemprate(j,h))
enddo
enddo
enddo
minwclose(2)=sum0/sum1
unempall(2)=(sum(px(8:n_x))-sum1)/sum(px(8:n_x))

!!!!! x specific bite !!!!!
do h=1,n_x
sum0=0.d0;sum1=0.d0
do i=1,n_a
do j=1,n_d

if (wage(j,h)<minw(h)*1.10d0) then
sum0=sum0+distw(i,j,h)*(1.-unemprate(j,h))
endif
sum1=sum1+distw(i,j,h)*(1.-unemprate(j,h))
enddo
enddo
minwclosexspec(h)=sum0/sum1
unempxspec(h)=(px(h)-sum1)/px(h)

enddo


!sum0=0.d0
!do jj=1,n_z
!sum0=sum0+(log(wage(jj,h))-log(wage(j,h)))/(log(zvect(jj))-log(zvect(j)))
!wageelast(j,h)=0.

!!!! Changing firms !!!!!
do i=1,n_a 
do j=1,n_d 
do h=1,n_x
ii=iprmat(i,j,h,2)
!temp=1.-unemprate(j,h)
sum0=0.d0

do jj=1,n_d
!sum0=sum0+px(h,hh)*pz(j,jj)*(1.-unemprate(jj,hh))*(pselect*choiceprobmat(ii,jj,hh)+(1.-pselect))
sum0=sum0+pselect*choiceprobmat(ii,jj,h)*(1.-pd(j,jj))*(1.-unemprate(jj,h))
enddo
changefirmrate(i,j,h)=sum0
enddo
enddo
enddo



!!!! Changing deciles !!!!
do i=1,n_a 
do j=1,n_d
do h=1,n_x
ii=iprmat(i,j,h,2)
jz=zcomp(j)
!temp=1.-unemprate(j,h)
sum0=0.d0
do jj=1,n_d
jjz=zcomp(jj)
if (decilez(jz).ne.decilez(jjz)) then
sum0=sum0+(pselect*choiceprobmat(ii,jj,h)+(1.-pselect)*pd(j,jj))
endif
enddo
changedecrate(i,j,h)=sum0

enddo
enddo
enddo

sum0=0.d0;sum1=0.d0;sum2=0.d0;sum3=0.d0;sum4=0.d0;sum5=0.d0
do i=1,n_a 
do j=1,n_d 
do h=1,n_x
ii=iprmat(i,j,h,2)
do jj=1,n_d
sum0=sum0+distw(i,j,h)*(1.-unemprate(j,h))*(1.-unemprate(jj,h))*&
     (pselect*choiceprobmat(ii,jj,h)+(1.-pselect)*pd(j,jj))*log(wage(j,h))*log(wage(jj,h))
     
sum1=sum1+distw(i,j,h)*(1.-unemprate(j,h))*(1.-unemprate(jj,h))*&
     (pselect*choiceprobmat(ii,jj,h)+(1.-pselect)*pd(j,jj))*log(wage(j,h))   

sum2=sum2+distw(i,j,h)*(1.-unemprate(j,h))*(1.-unemprate(jj,h))*&
     (pselect*choiceprobmat(ii,jj,h)+(1.-pselect)*pd(j,jj))*log(wage(jj,h))

sum3=sum3+distw(i,j,h)*(1.-unemprate(j,h))*(1.-unemprate(jj,h))*&
     (pselect*choiceprobmat(ii,jj,h)+(1.-pselect)*pd(j,jj))*log(wage(j,h))**2.d0   

sum4=sum4+distw(i,j,h)*(1.-unemprate(j,h))*(1.-unemprate(jj,h))*&
     (pselect*choiceprobmat(ii,jj,h)+(1.-pselect)*pd(j,jj))*log(wage(jj,h))**2.d0  

sum5=sum5+distw(i,j,h)*(1.-unemprate(j,h))*(1.-unemprate(jj,h))*&
     (pselect*choiceprobmat(ii,jj,h)+(1.-pselect)*pd(j,jj))    
enddo


enddo
enddo
enddo


autocovw=sum0/sum5-(sum1/sum5)*(sum2/sum5)
autocorrw=autocovw/((sum3/sum5-(sum1/sum5)**2.d0)*(sum4/sum5-(sum2/sum5)**2.d0))**0.5d0

!write(*,*) autocovw,autocorrw,sum1,sum2,sum3,sum4,sum5


sum0=0.d0;sum1=0.d0;sum2=0.d0;sum3=0.d0;sum4=0.d0;sum5=0.d0
sum6=0.d0;sum7=0.d0;sum8=0.d0;sum9=0.d0;sum10=0.d0
do j=1,n_d 
do jj=1,n_d
jz=zcomp(j)
jtau=taucomp(j)

sum0=sum0+distfirm(jz,jtau)*pd(j,jj)*log(yfirm(j))*log(yfirm(jj))    
sum1=sum1+distfirm(jz,jtau)*pd(j,jj)*log(yfirm(j))   
sum2=sum2+distfirm(jz,jtau)*pd(j,jj)*log(yfirm(jj))
sum3=sum3+distfirm(jz,jtau)*pd(j,jj)*log(yfirm(j))**2.d0   
sum4=sum4+distfirm(jz,jtau)*pd(j,jj)*log(yfirm(jj))**2.d0  

sum5=sum5+distfirm(jz,jtau)*pd(j,jj) 


sum6=sum6+distfirm(jz,jtau)*pd(j,jj)*log(empfirm(j))*log(empfirm(jj))    
sum7=sum7+distfirm(jz,jtau)*pd(j,jj)*log(empfirm(j))   
sum8=sum8+distfirm(jz,jtau)*pd(j,jj)*log(empfirm(jj))
sum9=sum9+distfirm(jz,jtau)*pd(j,jj)*log(empfirm(j))**2.d0   
sum10=sum10+distfirm(jz,jtau)*pd(j,jj)*log(empfirm(jj))**2.d0  
 
enddo
enddo


autocovy=sum0/sum5-(sum1/sum5)*(sum2/sum5)
autocorry=autocovy/((sum3/sum5-(sum1/sum5)**2.d0)*(sum4/sum5-(sum2/sum5)**2.d0))**0.5d0

autocovemp=sum6/sum5-(sum7/sum5)*(sum8/sum5)
autocorremp=autocovemp/((sum9/sum5-(sum7/sum5)**2.d0)*(sum10/sum5-(sum8/sum5)**2.d0))**0.5d0




open(unit=905,file="labor.txt")
    do j=1,n_d
    write(905,*) labor(j,1),labor(j,2)
    enddo
close(905)


open(unit=907,file="pi.txt")
    write(907,*) pitot
close(907)

open(unit=908,file="distw.txt")
    do i=1,n_a
    do j=1,n_d
    write(908,*) (distw(i,j,h),h=1,n_x)
    enddo
    enddo
close(908)

open(unit=909,file="wage.txt")
open(unit=910,file="unemp.txt")
    do j=1,n_d
    write(909,*) (wage(j,h),h=1,n_x)
        write(910,*) (unemprate(j,h),h=1,n_x)
    enddo
close(909)
close(910)

open(unit=911,file="f1.txt")
    do i=1,n_a
    do j=1,n_d
    write(911,*) (f1(i,j,h),h=1,n_x)
    enddo
    enddo
close(911)

open(unit=9112,file="xvect.txt")
    do h=1,n_x
    write(9112,*) xvect(h)
    enddo
close(9112)

open(unit=903,file="load_theta.txt")
    write(903,*) theta
close(903)

open(unit=292,file="pol.txt")
do i=1,n_a
do h=1,n_x
do j=1,n_d
jz=zcomp(j);jtau=taucomp(j)
!write (292,'(F18.8,A,F15.8,A,F15.8,A,F15.8,A,F18.10,A,F15.7,A,&
!&F20.8,A,F20.8,A,F20.5,A,F12.5,A,F15.8,A,F15.13,A,F15.13,A,&
!&F15.10,A,F18.8,A,F12.10,A,F12.10,A,F12.10,A,F20.10,A,F20.10,A,F17.10,A,I2,A,F20.8,A,F20.8)') &
!avect(i),',',zvect(jz),',',tauvect(jtau),',',xvect(h),',',wage(j,h),',',unemprate(j,h),',',&
!labor(j,1),',',laborfirm(j,1),',',kfirm(j),',',yfirm(j),',',empfirm(j),',',distw(i,j,h),',',distfirm(jz,jtau),',',&
!choiceprobmat(i,j,h),',',avect(iprmat(i,j,h,2)),',',changedecrate(i,j,h),',',mudem(j,h)&
!,',',changefirmrate(i,j,h),',',vsum(i,h),',',vinit(i,j,h),',',minw(h),',',h,',',&

write (292,'(F18.8,A,F15.8,A,F15.8,A,F15.8,A,F18.10,A,F15.7,A,&
&F20.8,A,F20.8,A,F20.5,A,F12.5,A,F15.8,A,F15.13,A,F15.13,A,&
&F15.10,A,F18.8,A,F12.10,A,F12.10,A,F12.10,A,F20.10,A,F20.10,A,F17.10,A,I2)') &
avect(i),',',zvect(jz),',',tauvect(jtau),',',xvect(h),',',wage(j,h),',',unemprate(j,h),',',&
labor(j,1),',',laborfirm(j,1),',',kfirm(j),',',yfirm(j),',',empfirm(j),',',distw(i,j,h),',',distfirm(jz,jtau),',',&
choiceprobmat(i,j,h),',',avect(iprmat(i,j,h,2)),',',changedecrate(i,j,h),',',mudem(j,h),',',changefirmrate(i,j,h),',',&
vsum(i,h),',',vinit(i,j,h),',',minw(h),',',h
enddo
enddo
enddo
close(292)

!open(unit=298,file="pol_choices_z.txt")
!do i=1,n_a
!do h=1,n_x
!do j=1,n_d
!do k=1,2
!jz=zcomp(j);jtau=taucomp(j)
!write (298,'(F18.8,A,F15.8,A,F15.8,A,F15.8,A,I2,A,F15.10,A,F15.10,A,F15.8,A,F15.8,A,F18.10,A,F15.7)') &
!avect(i),',',zvect(jz),',',tauvect(jtau),',',xvect(h),',',k,',',&
!choiceprobmat(i,j,h),',',avect(iprmat(i,j,h,k)),',',vsum(i,h),',',vlast(i,j,h,k),',',wage(j,h),',',unemprate(j,h)
!enddo
!enddo
!enddo
!enddo
!close(298)
!
!
!open(unit=298,file="pol_choices.txt")
!do ii=1,n_a
!do i=1,n_a,3
!do h=1,n_x,5
!do j=1,n_d,3
!do k=1,2
!jz=zcomp(j);jtau=taucomp(j)
!write (298,'(F18.8,A,F15.8,A,F15.8,A,F15.8,A,F15.8,A,I2,A,F15.10,A,F15.8,A,F15.8,A,F15.8)') &
!avect(ii),',',avect(i),',',zvect(jz),',',tauvect(jtau),',',xvect(h),',',k,',',&
!choiceprobmat(i,j,h),',',avect(iprmat(i,j,h,k)),',',vfun(ii,i,j,h,k),',',vlast(i,j,h,k)
!enddo
!enddo
!enddo
!enddo
!enddo
!close(298)


open(unit=209,file="keqnfinal.txt")
write (209,*) ksup,kdem,r
close(209)

open(unit=29,file="moments_fortran.txt")
!write (29,*) sdyfirm
!write (29,*) sdempfirm
!write (29,*) sdwage
write (29,*) unempall(1),',',unempall(2)
write (29,*) minwclose(1),',',minwclose(2)
write (29,*) autocorry
write (29,*) autocorremp
write (29,*) autocorrw
write (29,*) minwavgwagecollar(1),',',minwavgwagecollar(2)
close(29)

open(unit=64,file="minwagexspec.txt")
    do h=1,n_x
    write(64,*) minwclosexspec(h),unempxspec(h),avgwage(h),minw(h)
    enddo
close(64)

open(unit=904,file="avgl.txt")
    do j=1,n_d
    write(904,*) avgl(j)
    enddo
close(904)

open(unit=904,file="empfirm.txt")
    do j=1,n_d
    write(904,*) empfirm(j)
    enddo
close(904)


!diffminw(1)=minwclose(1)-0.05d0
!diffminw(2)=minwclose(2)-0.043d0

diffminw(1)=0.5d0*(minwavgwagecollar(1)-0.66d0)
diffminw(2)=0.5d0*(minwavgwagecollar(2)-0.50d0)


if (iterminw<n_iterminw.and.sum(abs(diffminw))>0.002d0) then

if (abs(diffminw(1))>0.002d0) minwvect(1)=minwvect(1)*(1-1.d0*diffminw(1))
if (abs(diffminw(2))>0.002d0) minwvect(2)=minwvect(2)*(1-1.d0*diffminw(2))

minwvect(1)=max(0.d0,minwvect(1))
minwvect(2)=max(0.d0,minwvect(2))


minw(1:7)=minwvect(1)
minw(8:n_x)=minwvect(2)



iterminw=iterminw+1

write(*,*) 'iterminw',iterminw,diffminw,minwvect(1),minwvect(2)
else

iterminw=100
endif

open(unit=90,file="minw.txt")
     write(90,*) minwvect(1)
     write(90,*) minwvect(2)
     write(90,*) minwvect(3)
     write(90,*) minwvect(4)
     write(90,*) minwvect(5)
close(90)

enddo !iterminw!



!!!!! SIMULATIONS !!!!!!!!
!!! We don't use simulations for alternative models. !!!!
!
!  
!      open (unit=501,file = 'randomz_10m.txt')
!      open (unit=502,file = 'randomx_10m.txt')
!      open (unit=503,file = 'randomk_10m.txt')
!      open (unit=504,file = 'randomchoose_10m.txt') 
!      open (unit=505,file = 'randomu_10m.txt') 
!      do i=1,n_id
!      do t=1,n_time-1
!      read(501,'(F15.14)') randomz(i,t)
!      read(502,'(F15.14)') randomx(i,t)   
!      read(503,'(F15.14)') randomk(i,t)  
!      read(504,'(F15.14)') randomchoose(i,t)
!      read(505,'(F15.14)') randomu(i,t)
!      enddo
!      enddo
!      close(501)
!      close(502) 
!      close(503) 
!      close(504) 
!      close(505)
!
!do id=1,n_id  
!
!!!! Variables in simulations !!!
!
!
!a_sim(id,1)=1
!z_sim(id,1)=int(n_z/2)
!!x_sim(id,1)=int(n_x/2)
!choose_sim(id,1)=1
!tau_sim(id,1)=int(n_tau/2)
!
!d_sim(id,1)=firmconvert(z_sim(id,1),tau_sim(id,1))
!
!x_sim(id,1)=1
!do while (sum(px(1:x_sim(id,1)))<randomx(id,1).and.x_sim(id,1)<n_x)
!x_sim(id,1)=x_sim(id,1)+1
!enddo
!
!enddo
!
!do id=1,n_id
!
!!print *, id
!if (id-int(id/100)*100==0) then
!print *, id
!endif
!
!do t=1,n_time-1
!if (choose_sim(id,t)==1) then !11
!
!k=1
!do while (sum(choiceprobmat(a_sim(id,t),1:k,x_sim(id,t)))<randomk(id,t).and.k<n_d)
!k=k+1
!enddo
!d_sim(id,t)=k
!z_sim(id,t)=zcomp(k);tau_sim(id,t)=taucomp(k)
!
!
!
!val_sim(id,t)=vsum(a_sim(id,t),x_sim(id,t))
!else !11
!!z_sim(id,t)=z_sim(id,t)
!val_sim(id,t)=vinit(a_sim(id,t),d_sim(id,t),x_sim(id,t))
!
!endif !11
!
!
!u_sim(id,t)=2
!if (randomu(id,t)<unemprate(d_sim(id,t),x_sim(id,t))) u_sim(id,t)=1
!
!
!if (u_sim(id,t)==2) then
!inc_sim(id,t)=wage(d_sim(id,t),x_sim(id,t))
!else
!inc_sim(id,t)=bmat(d_sim(id,t),x_sim(id,t))
!endif
!
!a_sim(id,t+1)=iprmat(a_sim(id,t),d_sim(id,t),x_sim(id,t),u_sim(id,t))
!
!c_sim(id,t)=inc_sim(id,t)+avect(a_sim(id,t))*(1.+r)-avect(a_sim(id,t+1))+pitot
!
!d_sim(id,t+1)=1
!do while (sum(pd(d_sim(id,t),1:d_sim(id,t+1)))<randomz(id,t).and.d_sim(id,t+1)<n_d)
!d_sim(id,t+1)=d_sim(id,t+1)+1
!enddo
!tau_sim(id,t+1)=taucomp(d_sim(id,t+1))
!z_sim(id,t+1)=zcomp(d_sim(id,t+1))
!
!!tau_sim(id,t+1)=tau_sim(id,t)
!!d_sim(id,t+1)=firmconvert(z_sim(id,t+1),tau_sim(id,t+1))
!
!x_sim(id,t+1)=x_sim(id,t)
!!x_sim(id,t+1)=1
!!do while (sum(px(x_sim(id,t),1:x_sim(id,t+1)))<randomx(id,t).and.x_sim(id,t+1)<n_x)
!!x_sim(id,t+1)=x_sim(id,t+1)+1
!!enddo
!
!choose_sim(id,t+1)=1
!if (pselect<randomchoose(id,t)) choose_sim(id,t+1)=0
!
!
!enddo !t!
!enddo !id!
!
!
!
!      open(unit=203,file="sim.txt")
!      do id=1,n_id
!      do t=1,n_time-1
!      write (203,'(I4,A,I4,A,F18.6,A,F18.6,A,F15.6,A,F15.6,A,I1,A,I1&
!&,A,F20.6,A,F20.6,A,F20.6,A,F18.6,A,F20.6,A,F20.6,A,F20.6,A,F20.6,A,F18.6,A,F15.6,A,I2)') &
!id,',',t,',',avect(a_sim(id,t)),',',zvect(z_sim(id,t)),',',tauvect(tau_sim(id,t)),',',xvect(x_sim(id,t)),',',u_sim(id,t),',',&
!choose_sim(id,t),',',val_sim(id,t),',',inc_sim(id,t),',',c_sim(id,t),',',minw(x_sim(id,t)),',',&
!laborfirm(d_sim(id,t),1),',',laborfirm(d_sim(id,t),2),',',kfirm(d_sim(id,t)),',',yfirm(d_sim(id,t)),',',empfirm(d_sim(id,t)),',',&
!distfirm(z_sim(id,t),tau_sim(id,t)),',',x_sim(id,t)
!     enddo  
!     enddo
!     close(203)

end program




